|
|
|
|
|
|
|
|
BVP Signal Analysis - A Complete Tour |
| Tags | pre-process☁bvp☁filter |
Physiological signals are informational resources containing a great potential of knowledge to be extracted, however, this knowledge is not directly accessible after signal acquisition.
In order to convert information into knowledge, physiological signals should be pre-processed, ensuring noise attenuation from unwanted sources and removal of unexpected artifacts. After pre-processing the signal, researcher can be totally secure that his valuable informational resource is ready to be processed, through the extraction of parameters.
Each one of the extracted parameters define an objective way to analyze the acquired and pre-processed physiological data. Subsequent interpretation of the results/extracted parameters will be extremely important to understand which type of events/mechanism are taking place inside the human organism and, at this point, an amazing journey is reaching its fundamental stage with the horizon of achievable discoveries becoming observable by the eyes of knowledge.
| ☌ The current Jupyter Notebook belongs to a set of Notebooks entitled "Complete Tour", where it is presented, in a sequential way, a guide containing recommended steps to be followed, namely: 1) sensor/electrode placement; 2) pre-acquisition notes; 3) pre-processing methods (filtering and motion artifact reduction); 4) detection of specific events and 5) parameter extraction procedures. |
A - Blood Volume Pulse (BVP) sensor | Finger Probe Placement
Taking into consideration the great density of capillaries and great translucency in the hand fingers, these become the ideal site to place the BVP sensor and our finger probe, as demonstrated in the illustrative example of the following figure:
B - Pre-Acquisition Requirements
In research literature it is considered that the typical informational band of photoplethismographic BVP signal is in the range:
-
[0.02; 2.1] Hz |
source link
- [0.6875; 10] Hz | source link
-
[0.5; 4] Hz |
source link
Accordingly to these articles, informational content of BVP signal reach, at the most broader criteria, the 10 Hz frequency components. So, applying the widely known
Nyquist Theorem
, before starting the data collection, the signal acquisition device must be configured to support sampling rates of at least
20 Hz
, avoiding aliasing problems during the analog-to-digital conversion procedure.
C - Pre-Processing Stage
C0 - Importing Packages and load signal
# biosignalsnotebooks own package.
import biosignalsnotebooks as bsnb
# Scientific programming package.
from numpy import average, std, diff, max, pad, log2, ceil, absolute, power
from numbers import Number
# Pacakge dedicated to Wavelet decomposition algorithms.
from pywt import swt, iswt
# Deep copy function to clone lists.
from copy import deepcopy
# Load entire acquisition data.
data, header = bsnb.load("../../signal_samples/bvp_motion_artifact.txt", get_header=True)
# Store the desired channel (CH2) in the "signal" variable
signal = data["CH2"]
# Sampling rate definition.
sr = header["sampling rate"]
C1 - Filtering Stage 1
In this filtering stage it will be used a conventional bandpass filter, intended to attenuate all signal components outside the typical BVP frequency bands.To compare the influence of multiple filter bandpass frequencies it will be analyzed a restrictive range - [0.02; 2.1] Hz - and a more extended option - [0.6875; 10] Hz.
C1.1 - Generation of filtered signals
# Conventional 2nd order Butterworth bandpass filtering
# [Option RES]
signal_res = bsnb.bandpass(signal, 0.02, 2.1, order=2, fs=sr, use_filtfilt=True)
# Conventional 2nd order Butterworth bandpass filtering
# [Option EXT]
signal_ext = bsnb.bandpass(signal, 0.6875, 10, order=2, fs=sr, use_filtfilt=True)
C1.2 - Comparison between filtered and unfiltered signals
With this first filtering stage we achieved a reasonable result:
- The restrictive filtered signal ( signal_res ) generates an excessive smoothed signal, causing lost of information, namely the dicrotic notch ;
- The extended filtered signal ( signal_ext ) was successful in smoothing the signal but the highlighted artifacts and hand movement oscillations were kept;
Taking this information into consideration, it will be necessary to include an additional filtering stage. As a consequence of the preservation of all relevant information, in the subsequent steps only the signal_ext will be analyzed.
The following section is dedicated to minimize the artifacts influence on our BVP signal.
C2 - Filtering Stage 2
After smoothing the BVP signal, removing high-frequency noise components, in the following steps it will be adapted an algorithm proposed by C. M. Lee et al. ( research article
) to minimize the influence of motion artifacts in the BVP signal.
C2.1 - Decomposition of signal_ext using Stationary Wavelet Transform (SWT)
SWT has the advantage, when comparing with Discrete Wavelet Transform (DWT), of not producing additional "ringing effects in the vicinity of a discontinuity", accordingly to Chen et al. ( research article
).
# SWT 7th level Decomposition using "Daubechies" mother wavelet [https://ieeexplore.ieee.org/abstract/document/6908199].
signal_ext = pad(signal_ext, (0, 2**int(ceil(log2(len(signal_ext)))) - len(signal_ext)), "constant") # Ensure that the number of signal samples is a power of 2.
swt_orig_coeffs = swt(signal_ext[:32768], "db4", level=7)
C2.2 - Removal of motion artifact wavelet coefficients (SWT scaled/approximation coefficients) between decomposition levels 1 and 7 (list entries 0 to 6)
A reference threshold (for detail coefficients) and a range (for scaling/approximation coefficients) are determined for each decomposition level, consisting in:
- the average value of the $N$ detail coefficients in the level under analysis (for detail coefficients $d_{2^j}^{2^jk + p}$, being $j$ the wavelet scale, $k$ the translation and $p$ the respective shift) \begin{equation} \mu_j^d = \frac{\sum_{n=1}^Nd_{2^j}^{2^jk + p}[n]}{N} \end{equation}
- a range delimited by lower and upper thresholds: $T_{low} = \mu_j^a - \alpha \times \sigma_j^a$ and $T_{high} = \mu_j^a + \alpha \times \sigma_j^a$, where $\alpha$ defines a positive constant factor, that controls the range width (bigger values will generate an extended range), $\mu_j^a$ is the average value of $N$ approximation coefficients ($a_{2^j}^{2^jk + p}$, being $k$ the translation and $p$ the respective shift) in the level $j$ under analysis while $\sigma_j^a$ is the respective standard deviation: \begin{equation} \mu_j^a = \frac{\sum_{n=1}^Na_{2^j}^{2^jk + p}[n]}{N}\quad\quad\sigma_j^a = \sqrt{\frac{\sum_{n=1}^N(a_{2^j}^{2^jk + p}[n] - \mu_j^a)^2}{N}} \end{equation}
The removed coefficients are related with motion artifacts, taking into consideration that they are outliers outside the acceptability ranges ($[\mu_j^d; +\infty]$ and $[T_{low}; T_{high}]$).
# "level" iteration for loop.
swt_coeffs_copy = deepcopy(swt_orig_coeffs)
for level in range(0, 7):
thr_avg_dt = average(swt_orig_coeffs[level][1]) # Only "detail" coefficients are being taken into account.
thr_avg_sc_low = average(swt_orig_coeffs[level][0]) - 3 * std(swt_orig_coeffs[level][0])
thr_avg_sc_high = average(swt_orig_coeffs[level][0]) + 3 * std(swt_orig_coeffs[level][0])
# "coefficients" iteration loop.
for coeff_nbr in range(0, len(swt_orig_coeffs[level][1])):
if swt_orig_coeffs[level][1][coeff_nbr] > thr_avg_dt: # Motion artifact coefficients.
swt_orig_coeffs[level][1][coeff_nbr] = 0
if swt_orig_coeffs[level][0][coeff_nbr] < thr_avg_sc_low or swt_orig_coeffs[level][0][coeff_nbr] > thr_avg_sc_high: # Motion artifact coefficients.
swt_orig_coeffs[level][0][coeff_nbr] = 0
else: # Storage of noise artifact coefficients in a separate list.
swt_coeffs_copy[level][0][coeff_nbr] = 0
C2.3 - Signal reconstruction through an inverse SWT
rec_signal = iswt(swt_orig_coeffs, "db4")
In spite of a considerable improve in the quality of our signal, a parcel of the artifact still exists.
Fortunately, the wavelet coefficients of the noise artifact was stored on step C2.2 in the swt_coeffs_copy variable.
Now we can reconstruct the noise artifact component.
C3 - Noise artifact isolation
C3.1 - Signal reconstruction
art_signal = iswt(swt_coeffs_copy, "db4")
C3.2 - Enhancement of noise complex
In order to isolate the noise artifact complex the following steps will be applied:- Signal rectification
- Signal smoothing with a big size moving average window (size equal to 2 times the number of samples acquired per second - sampling rate), because we want an envelope instead of the details
- Increase of the separation between variant and quasi-invariant components of the signal (generating the 2nd power signal - this order can be customized)
# Signal rectification.
rect_art_signal = absolute(art_signal)
# Signal smoothing stage.
smooth_art_signal = bsnb.smooth(rect_art_signal, 2 * sr)
# Increase separation between variant and quasi-invariant components of the signal.
power_art_signal = power(smooth_art_signal, 2)
C3.3 - Signal binarisation
All sample values greater than 1 (this value may not be appropriate for all BVP acquisitions) will be equal to 1 while the remaining samples will be updated to 0.noise_binary_complex = []
for i in range(0, len(power_art_signal)):
if power_art_signal[i] > 1:
noise_binary_complex.append(1)
else:
noise_binary_complex.append(0)
C3.4 - Identification of the start and end of each noise artifact complex
Upward transitions define the start of noise complexes while downward transitions are the respective ends.diff_signal = diff(noise_binary_complex)
noise_complexes = [(None, None)]
for i in range(0, len(diff_signal)):
# Search for a starting period.
if diff_signal[i] == 1 and noise_complexes[-1][0] == None:
# Convert tuple to list.
tuple_list = list(noise_complexes[-1])
# Change entry reserved to the "end" value.
tuple_list[0] = i
# Update noise_complexes content.
noise_complexes[-1] = tuple(tuple_list)
# Search for an ending period.
if diff_signal[i] == -1 and noise_complexes[-1][1] == None and noise_complexes[-1][0] != None:
# Convert tuple to list.
tuple_list = list(noise_complexes[-1])
# Change entry reserved to the "end" value.
tuple_list[1] = i
# Update noise_complexes content.
noise_complexes[-1] = tuple(tuple_list)
# Preparing for the next searching iteration.
noise_complexes.append((None, None))
# Remove last tuple if it does not have an end.
if noise_complexes[-1][1] == None:
noise_complexes.pop()
C4 - Removal of the remaining influence of noise complexes from the BVP signal
for i in range(0, len(noise_complexes)):
noise_period_start = noise_complexes[i][0]
noise_period_end = noise_complexes[i][1]
rec_signal[noise_period_start:noise_period_end] = average(rec_signal)
D - Parameter Extraction
From the filtered BVP signals it is possible to execute a deep heart rate variability analysis, similar to the one presented in another Jupyter Notebook ( ECG Analysis - Heart Rate Variability Parameters"
).
A more specific analysis focused on the BVP morphology will be soon available, but, for now, we invite you to explore a very interesting research article from Mohamed Elgendi (
source link
).
With the steps described on the current Jupyter Notebook you are able to acquire a broad perspective of the entire methodology involved on BVP signal analysis, starting in the acquisition stage and reaching the parameter extraction.
This is a mere illustrative example aiming to show some possible paths, but, you are free to explore this amazing world by your own.
We hope that you have enjoyed this guide.
biosignalsnotebooks
is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining
Notebooks
!